// Copyrights Orestis Georgiou   25-8-09
//
// This program is related to the paper: 
//Product of n independent Uniform Random Variables.

// When Theorem 1. of the above paper is integrated, the incomplete gamma 
// function, Gamma[n,x] is obtained which for integer n can be expressed
// analytically. Thus, this program calculates the probability that a random 
// variable X, consisting of the product of n independent and identically
// distributed uniform [a,b] random variables X_{i}, i=1,2..n, takes on a value
// less than or equal to tau.
//
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// define the endpoints of the supports of each function
#define alpha(k) (pow(a,(n-k+1))* pow(b,(k-1)))
#define beta(k) ( pow(a,(n-k))* pow(b,k))

 int n=4;  // input here the value of n
 
 int j, m;
 double a= 0.3, b=3, tau= 0.31; // input here the values of a, b and tau
 double Prob, diff, OutPut, A;
// define the factorial function
double fac( int n )
{double fact = 1.0;   
   while ( n > 1) {
       fact = fact * n;
       n = n - 1;}
   return fact;}
// define the F2 function
double F2(int k)
{
for(j=0,OutPut=0;j<(n-k+1);j++)
   {
   for(m=1,A=0;m<n;m++)
             {A+= pow(log( (pow(a,j)*pow(b,n-j))/beta(k) ), m) / fac(m);}
   OutPut+=beta(k)*fac(n-1)*(1+A)*n*pow((-1),j)/(fac(j)*fac(n-j)*pow((b-a),n));
   }
return OutPut;
}
// define the F1 function
double F1(int k)
{
for(j=0,OutPut=0;j<n-k+1;j++)
   {
   for(m=1,A=0;m<n;m++)
             {A+= pow(log( (pow(a,j)*pow(b,n-j))/alpha(k) ), m) / fac(m);}
   OutPut+=alpha(k)*fac(n-1)*(1+A)*n*pow((-1),j)/(fac(j)*fac(n-j)*pow((b-a),n));
   }
return OutPut;
}
// define the Ftau function
double Ftau(int k)
{
for(j=0,OutPut=0;j<n-k+1;j++)
   {
   for(m=1,A=0;m<n;m++)
             {A+= pow(log( (pow(a,j)*pow(b,n-j))/tau ), m) / fac(m);}
   OutPut+= tau*fac(n-1)*(1+A)*n* pow((-1),j)/(fac(j)*fac(n-j)*pow((b-a),n));
   }
return OutPut;
}
// the main code!!!
int main()
{
  int k=1;
  
  Prob = 0.00;
  while ( tau > beta(k)) {            
        diff= F2(k)-F1(k);
        Prob= Prob + diff;
        k=k+1;       
   }
Prob+= Ftau(k) - F1(k);
printf("\na= %g,\t b= %g,\t tau= %g,\n\nP(X <= tau) = %.17f \n\n",a,b,tau,Prob);
system ("PAUSE"); // remove this line according to your compiler
} 
//The END